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Dynamical Second Order Closure Theories 


2.1 Stochastic averaging, non-equilibrium 
statistical mechanics, and quasilinear 
approaches — Bouchet et al. 

2.2 Stochastic structural stability theory — 
Farrell & Ioannou 

2.3 Zonostrophic instability — Young &; 
Srinivasan 


2.4 Zonal Flow as Pattern Formation 
Parker and J. A. Krommes 


J. B. 


2.4.1 Introduction 

This section continues the use of statistical methods to in¬ 
vestigate the physics of zonal jets. Our interest in the prob¬ 
lem of turbulent-driven zonal flows stems from wanting to 
understand their behavior in magnetized plasmas. Plasmas 
possess their own host of complexities distinct from those of 
geophysics, including the mass differences of ions and elec¬ 
trons, kinetic effects such as wave-particle resonances, and 
electromagnetic effects. It is a marvel that despite the im¬ 
mense disparities between laboratory plasmas and planetary 
atmospheres, similar physics in each conspire to organize 
regular flows out of turbulence. 

In the challenging envir onments of plas mas devices, zonal 
flows have been detected ( Fuiisawal . 21 It 1 9 i. Likewise, sophis¬ 
ticated gyrokinetic simulations, which are thought to retain 
all of the important physics of plas ma microt urbulence, ex¬ 
hibit the formation of zonal flows ( Lin et all Il998l b Zonal 
flows have taken on special significance in magnetized plas¬ 
mas because these flows are thought to suppress turbulent 
transport of heat ( Diamond et all I2005T ). an advantageous 
feature when the ultimate goal is to keep hot the core plasma 
of fusion reactors. As a result, much research into the physics 
of zonal flows in plasmas has been undertaken, and a great 
deal has been learned (see sections 3.2.1 and 3.2.2). 

But several questions remain. One problem of basic inter¬ 
est involves the length scale of the zonal flows (the width 
of the jets). No one has yet found a heuristic estimate of 
the jet width in plasmas that enjoys as much success as the 
Rhines scale in geophysical contexts. Another question is the 
detailed mechanism by which the flows suppress turbulence 
( Hatch et al. I [20131) . The flows can stabilize the linear modes 
responsible for driving the turbulence. Additionally, plasma 


possess a complex array of feedbacks as well as velocity-space 
structure, any of which might be responsible for dissipation. 
Strong interactions between turbulence and zonal flows can 
modify the energetics. This is an active research area which 
requires strong collaboration between experiment, computa¬ 
tion, and theory. 

With that in mind, using the simplest paradigm systems 
we seek to develop a basic theory of zonal flows that can 
serve as the foundation for more complete plasma models. 
To that end we begin from the Hasegawa-Mima equation, 
which describes electrostatic plasma turbulence in a uniform 
magnetic fiel d in the presence of a backgr ou nd plasma den ¬ 
sity gradient IlH asegawa and Mimal. 1 19781 : ISmolvakov et all . 
hoQQal : iKrommes and Kiml.l200dj . This equation is equiva¬ 
lent to the quasigeostrophic barotropic vorticity equation in 
a certain limit. 

We have found that from a statistical perspective, zonal 
flows constitute pattern for mation amid a turbulent bath 
( Parker and Krommed . l2013bll . Our account here empha¬ 
sizes the role of a symmetry breaking and its consequences. 
Some of the key insights to emerge are a mathematical pre¬ 
diction of a nonunique jet wavelength and a close linking of 
the phenomenon of jet merging with stability of the zonal 
flow-turbulence equilibrium. We also study the symmetry¬ 
breaking zonostrophic instability in some detail and add 
some novel insights. 

We are mindful that the Hasegawa-Mima equation is not 
a realistic or quantitatively accurate model of magnetized 
plasmas. We are not primarily concerned with specific pa¬ 
rameter dependencies, but rather we wish to establish gen¬ 
eral principles that will act as the building blocks in more 
elaborate theories and models. The Hasegawa-Mima equa¬ 
tion is a minimal model in that it contains the necessary 
physics for zonal flows to form but other complicating de¬ 
tails are stripped away. 

The zonal flows that are generated in these simple models 
can in some regimes be steady in time, or at least nearing 
an idealization where that is true. Such steady zonal flows 
may in fac t occur i n nature, such as the je ts in Jupiter’s at¬ 
mosphere l Yasavada and Showmanl . hoosl b Steady jets may 
also be a robust fea ture of plasma turbulence in uniform 
magnetic geometry ( Numata et all l2007lb In the toroidal 
geometry relevant to fusion plasma, the zonal flows may 
fluctuate in time. Nevertheless, for a tractable starting point 
we restrict ourselves to consideration of zonal flows that are 
steady or perhaps evolving much more slowly than the tur¬ 
bulence. Although the Hasegawa-Mima equation in some 
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regimes produces nonsteady jets, we select a parameter 
regime where the jets are nearly steady. In any realistic situ¬ 
ation there will always be some variation in time, and finding 
a true steady state requires a statistical perspective. 

We have adopted the statistical approach to understand¬ 
ing turbulence. While a brief review is given here, a broad in¬ 
troduction to statistical turbulence is given in Section 5.1.1. 
This approach complements other methods such as making 
detailed measurements of plasma fluctuations or perform¬ 
ing direct numerical simulations (DNS), which accumulate 
reams of data so vast that it can be unclear how one should 
go about making sense of it all. The aim of the statisti¬ 
cal approach is to focus on the macroscopic quantities of 
interest, such as transport coefficients, energy spectra, and 
the like. By working with averaged quantities from the out¬ 
set, one can circumvent the rapid spatiotemporal fluctua¬ 
tions and potentially see a clearer view of the physics. Of 
course, there is no free lunch. As a consequence of averaging 
a nonlinear equation, one is generally left with the average 
of an unknown quantity: a closure problem. Various statisti¬ 
cal closures, perhaps motivated by physical considerations, 
provide different approximations for the unknown terms. A 
major difficulty with this approach is that the closures are 
essentially uncontrolled approximations; the nonlinearity in¬ 
herent to turbulence makes it hard to know exactly what 
is lost. The closure might obliterate some highly coherent 
or correlated phenomena. Nevertheless, these difficulties do 
not invalidate th e statistical approach, from which much 
has been le arned ( Frischl . fl995l : iKrommesl . 120021 : iKraichnanl . 
1959, 1964). Historically, the majority of theoretical studies 
into turbulence that follow this approach assume homoge¬ 
neous statistics, where the statistics of turbulent quantities 
do not depend on position. Consequently, most of the theo¬ 
retical machinery that has been developed also applies only 
to homogeneous statistics, with comparatively little devoted 
to inhomogeneous statistics. 

In the presence of steady zonal flows, one is inevitably led 
to the conclusion that a proper statistical description must 
allow for inhomogeneous turbulence. The turbulence cannot 
be statistically homogeneous; a location within the peak of a 
jet differs physically from a location at the node of a jet. In¬ 
homogeneous turbulence often arises on large spatial scales 
due to the presence of boundaries, topography, inhomoge¬ 
neous driving forces, or other external factors. These may 
be present even though it is often assumed that turbulence 
homogenizes on small scales. But external influence is not 
the only way that inhomogeneities might develop. Even if 
topography and boundaries are removed and the problem is 
contrived such that the governing equations of motion are 
independent of position, that translational symmetry may 
be broken spontaneously. Zonal flows and inhomogeneous 
turbulence can result from such a scenario. In this subtler 
development of inhomogeneous turbulence where the under¬ 
lying physics is translationally invariant, a satisfactory sta¬ 
tistical description must still allow for inhomogeneity. 

As with several other sections in this chap¬ 
ter, we e mploy the se co nd-or de r cumulant (CE2) 
framework dM arston e t all 20081 : Tobias et al. , 2011 


iTobias and Marstonl . 2014? iFarrell and Ioannoul . l2003l . 120071 . 


2009 : 

2013; 


^akas_andHoannou, 2011. 20 13bl : IConstantinou et all 
Parker and Krommej. l2013bli Csee also sections 5.1.1, 


5.1.2, 5.2.2, 5.2.3, and 5.2.4). The CE2 formalism is the 
simplest possible setting in which to study inhomogeneous 
turbulence with statistical equations. One way of deriv- 
in g CE2 is throug h t he qu asilinear (QL) approximation 
( Srinivasan and Yound . 120121 1. 

We remark that the CE2 framework is equivalent to the 
Wigner-Moyal formalism. The Wigner-Moyal formalism has 
been used in stu dies o f wave p h ysics in inhomogeneou s me¬ 
dia ( McDonald and Kaufmari 1 19851 : lHall et all l2002h . The 
Wigner distribution function, assuming an appropriate av¬ 
erage is used in its definition, is closely related to the CE2 
correlation function: they are both the two-point, one-time, 
second-order correlation of fluctuations. The Wigner-Moyal 
equation, which describes the evolution of the distribution 
function, is the analog of the CE2 equation. 

In this article we use the 2D Charney-Hasegawa-Mima 
equation (CHME), written in a form to also encompass the 
modified Hasegawa-Mima equation (mHME). The funda¬ 
mental equation i ^3 


dtw + v ■ Vic + fldxip = f + D 


( 2 . 1 ) 


where w is potential vorticity, v = z x Vi/’ is the hori¬ 
zontal velocity, ip is the stream function such that w = 
V 2 ip — azFLj 2 ^- The standard CHME is obtained with 
&zf = 1- The mHME involves setting &zf = 0 only for 
modes with k x = 0 , i.e., zonal flow modes, and &zf = 1 f° r 
all other modes. This is done to more accurately model the 
zonal flow response in plasmas. Additionally, in the plasma 
context, Ld is the sound gyroradius p s - f is the external forc¬ 
ing and can be thought of as a stirring, e.g., some idealiza¬ 
tion of excitation by buoyant convection or thermal gradient 
instability. D is dissipation and might represent frictional or 
viscous damping, or more generally any net transfer to exter¬ 
nal degrees of freedom. Since the system is both driven and 
damped, it reaches an equilibrium where energy injection is 
balanced by energy dissipation. 

To illustrate the QL approximation, we temporarily set 
/ and D to zero. Conceptually, one decomposes the flow 
field into a zonally symmetric part (the zonal flow) and the 
residual (the eddies or tubulence). Let w = w + w', where 
the overbar represents a zonal average, or average over x. 
Equation (ED can be decomposed as 


dtM1 + v' • Vw' = 0, (2.2a) 

dtw' + v ■ Vw 1 + v ■ Vw 

+ v' ■ Vw 1 — v' • Vw' + /3 dxtp' = 0. (2-2b) 

No approximation has been made thus far. At this point, one 
can make certain approximations that treat the eddies and 
the zonal flows differently. The QL approximation proceeds 
by dropping, within the eddy equation <lT2bl) . the terms 
quadratic in the eddy quantity (the advective nonlinearity). 

1 To obtain the coordinates conventionally used in the plasma 
literature, let {x,y} —> {—y,x}. 
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The QL system is 

d t w + v' • Vw' = 0, (2.3a) 

dtw' + v • Vw' + v / ■ Vu7 + /39x'0 , = 0. (2.3b) 

An alternative way of thinking about the QL approxima¬ 
tion is in Fourier space. All triad interactions between three 
Fourier modes are neglected except for those triads that con¬ 
tain one zonally symmetric mode ( k x = 0). The QL system 
respects the nonlinear conservation of energy and enstrophy. 
It should be noted, however, that the QL approximation de¬ 
stroys exact material conservation of potential vorticity. 

One of the uses of the QL system is that, subject only 
to an ergodic assumption, a statistical description can be 
obtained without a closure problem. Thus can one obtain 
CE2. In the CE2 framework, the dynamical variables are 
the zonal-mean field w and the two-point covariance W = 
(w'(xi)w'(x 2 )). The angle brackets may refer to a zonal 
average, ensemble average, or some other appropriate oper¬ 
ation. The derivation of CE2 from QL ensures that CE2 is 
statistically realizable with well-behaved statistics. 

If the flow is predominantly zonal, then the QL approx- 
i mation may b e valid, at least for eddies at large scales 
( Bouchet et, all hoiflh . However, we do not wish to restrict 
ourselves to discussion of a particular regime where zonal 
flow dominates. In fact, most of our work focuses on the 
opposite limit where the zonal flow is weak relative to the 
turbulent flow. We argue that even though the QL system 
is not rigorously valid as an approximation, it is useful as 
a model which contains some of the same behavior as the 
true system. In particular, numerical evidence shows that 
the same symmetry breaking occurs in the QL system as 
in the original system. The QL model is simpler and more 
tractable, however, and so provides a window into under¬ 
standing the physics. 

Although there are quantitative difference between CE2 
and the original nonlinear dynamical system, CE2 does have 
something useful to offer about the physics of zonal flows. 
CE2 provides a tractable problem with which to gain fun¬ 
damental insight into the behavior of zonal flows and their 
interaction with turbulence. 


2.4-1.1 Analogy between zonal flows and 
Rayleigh-Benard convection rolls 


The notion of spontaneous symmetry breaking with 
respect to zonal flows has been discussed befor e 
( Farrell and Ioannoul . 12007) : ISrinivasan and Younel l2012l l. 
This section will expand on that in discussing the mechanics 
of the symmetry breaking, as well a s specific consequences 
i t has for the physics of zonal flows (Parker and Krommes, 
l2013bllal ). 


An important aspect of zonostrophic instability is that 
it involves a spontaneous symmetry breaking. A sponta¬ 
neous symmetry breaking occurs when a situation’s govern¬ 
ing physics are invariant under a symmetry transformation 
but a physical realization is not invariant under the same 
transformation. A simple example would be a ball moving 
in a symmetric double-well potential, as in Figure l2Tl The 



Figure 2.1 Discrete spontaneous symmetry breaking occurs 
when a ball moving in a symmetric double-well potential must, 
due to friction, end up in one of the wells. 



Figure 2.2 Convection rolls in Rayleigh-Benard convection 
break the horizontal translational symmetry. 


equations of motion of the ball are invariant to reflection 
about the center line. But with friction the ball must even¬ 
tually end up in one of the wells, a state which breaks the 
symmetry. 

Another well-known example of spontaneous symmetry 
breaking is the for matio n of convection rolls in Rayleigh- 
Benard convection ( Busse . 1978 1. A box of fluid, taken to 
be infinite in both horizontal directions and finite in verti¬ 
cal extent, is heated from below. At weak heating, the heat 
is transferred to the cooler top surface solely by conduc¬ 
tion, and the fluid is motionless. But at sufficiently high 
heating, buoyancy forces overcome the inherent dissipation 
and the conduction state becomes unstable to the forma¬ 
tion of convection rolls, as shown schematically in Figure 
E3 The convection rolls are spatially periodic but steady 
in time. This transition to convection is analogous to the 
generation of zonal flows out of homogeneous turbulence. 
Like the conduction state, homogeneous turbulence is (sta¬ 
tistically) uniform in space. And as a drive parameter such 
as the strength of the forcing is varied, that uniform state 
becomes unstable to the formation of a periodic structure. 
Born out of turbulence are spatially periodic, steady-in-time 
zonal flows, which are analogous to the convection rolls (see 
Figure 12751) . More than merely descriptive, this analogy will 
be made mathematically precise in section 12.4.3. II 
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Figure 2.3 Zonal flows on a beta plane break the north-south 
(statistical) translational symmetry. 


2-4-1.2 Outline 

This rest of this section is separated into two main parts. 
The first part reexamines zonostrophic instability with the 
goal of improving physical understanding of the generation 
of zonal flows. We show that zonostrophic instability con¬ 
tains the 4-mode modulational instability as a special case. 
In |Appcndix 2.47A| we provide a physical picture of the in¬ 
stability in the limit of long-wavelength zonal flows. 

The second part studies the equilibrium between turbu¬ 
lence and zonal flows. The symmetry breaking and the bifur¬ 
cation to zonal flows is studied in some detail. In addition, 
one method of numerical solution to the CE2 equations is 
offered. 


2.4.2 Zonal Flow Generation through Instability of 
Homogeneous Turbulence 


ISrinivasan and Younel ( 2012 i ) gave a detailed and insightful 
treatment of the so-called zonostrophic instability. In zonos¬ 
trophic instability, a homogeneous turbulent background is 
unstable to coherent zonal flow perturbations. We give a 
brief overview before studying specific cases. Srinivasan and 
Young began with a convenient, real-space form of the CE2 
statistical equations. We generalize their work to allow for 
finite L^ and to unify the modified Hasegawa-Mima equa¬ 
tion and the (equivalent) barotropic vorticity equation. The 
appropriate CE2 equations, which can be derived from the 
QL equations (12.31) . are 


covariance of stream function and is given by 


L = 


(V + Oydy + ■ 


L't>{x,y | y,t), 

(2.5) 

<1 

to 

1 

+ 

1— 1 

(2.6) 


For simplicity the visc osity has been tak en to be zero ; though 
that is not necessary ( Parker and Krommesl . l2013bl a|) . The 
external forcing / has been taken to be random white noise, 
and F(x,y) is its covariance. 

One equilibrium of Eq. (i2~n> has no mean zonal flow, 
U = 0, and corresponds to homogeneous turbulence. The 
covariance is independent of y and takes the simple form 


W H (x,y) = 


F(x,y ) 

2p 


(2.7) 


This is always a steady-state solution of Eq. (l 2 Tl> . But that 
does not mean it will naturally occur; it may be unstable. 
Particularly of interest is the stability to perturbations with 
a mean-field component, i.e., zonal flow. These perturba¬ 
tions are assumed to have dependence such that A 

is the eigenvalue and q is th e wavenumber of t he zonal flow. 
The dispersion relation is ( Srinivasan and Yound . [20121 ) 


%(A + /x) =qA- — qA_|_, 

<r 


where q 2 = q 2 + a zfLj 2 , 


A± = 


/ 


dkxdky k 2 k y ( 1 - q 2 /h±)W H (h x , k v ± Iq) 
(2tt) 2 


(A + 2 y)h+h- + 2i/3qk x k y 


(2.8) 


(2.9) 


—2 q / 1 \2 o 

h.± = k x + (ky ± ^q) + L 7 , k x and k y are the Fourier 

conjugate variables of x and y , and our Fourier transform 
convention is 

f{k) = J dx e~ lkx f(x). ( 2 . 10 ) 

Although this convention uses the same symbol / for the 
real space and Fourier space functions, we always specify 
the argument of the function to make clear whether the real 
or Fourier domain is being used. In Eq. (12.81) . the left-hand- 
side (LHS) is the zonal flow intrinsic response and the right- 
hand-side (RHS) is the Reynolds stress forcing term. If the 
perturbations are unstable, corresponding to a solution with 
Re A > 0, then zonal flows grow in the so-called zonostrophic 
instability. One can also let p and F be zero; in that case, 
with U — 0 any homogeneous Wjj is a steady state and has 
the dispersion relation above. 


dtW + ( U+ - U-)d x W - (U+ - lf_) (V + i«9|) d x V 

- [2P - (U++ u'L)\dyd x d v q = F(x, y) - 2 gW, 

(2.4a) 

dtiu + dyd x d y <f{ 0 , 0 I y,t) = -nU, (2.4b) 

where W(x,y \ y,t) is the 2-point covariance of vorticity, x 
and y are difference coodinates and y is the average coordi- 

■ 2 Q _ 9 

nate, V = V — , U(y, t ) is the zonal-mean zonal veloc¬ 

ity, y is the scale-independent friction, U± = U(jj± ^y,t ), 
U± = U£-a ZF L- 2 U±,7 = l-azpL^dF 2 , and T is the 


2.4-2.1 Isotropic Background Spectrum with Finite 
Deformation Radius 


We now specialize the dispersion relation to an isotropic 
background and examine various limits, allowing for finite 
deformation radius. Although a purely isotropic spectrum is 
unlikely to obtain in practice when the beta effect is present, 
such an investigation helps to isolate the physical conse¬ 
quences of various effects. 

In the context of an infinite deformation radius L the 
effect of an isotropic backgroun d spect r um has been studied 
before ( Srinivasan and Younel 120121 : iBakas and loannoul . 
































l2013ah . Those studies concluded that for an isotropic back¬ 
ground, /3 7 ^ 0 is required for instability. Additionally, they 
found that for an isotropic background, the eddies acted on 
long-wavelength zonal flows as a negative hyperviscosity in¬ 
stead of negative viscosity. That is, the eddy forcing on the 
RHS of (12.81) behaves as g 4 rather than q 2 at small q. In 
this section, we study how these results change when finite 
deformation length Lj is allowed. 

The dispersion relation for a homogeneous, isotropic back¬ 
ground spectrum can be written as 

^(A + m) = i J ^ k 2 (l- S) W H (k)S( X ,n,m) 

( 2 . 11 ) 

where 


p2lT 

S{x,n,m) = / 

•to 


d(j) 


K, 


K = 


2n 

(n — 2 cos (j>) sin 2 (j> 


( 2 . 12 ) 


x(l — 2 n cos <j> + n 2 + m) + i(n — 2 cos rf>) sin <j> ’ 

(2.13) 

(A + 2 fijk 2 


X= to 

"=!■ 
m = ( kL d )~ 


(2.14) 

(2.15) 

(2.16) 


We now examine the limit of large X-, which could correspond 
to either small /3 or small q. Asymptotic expansion of S for 
large x reveals interesting behavior that can differ for finite 
vs. infinite Z,^. 

For infinite Lj (i.e., m = 0 ), S behaves as0 

( Srinivasan and Younel . I2OI2I ) 


S(x,n, 0) = < 


+ 0 (x“ 5 ) 


X 3 8(1 — n 2 ) 


1 n 2 - 1 


X 


2 n 3 


+ 0(x 3 ): 


n 2 < 1 , 


n z > 1 . 


(2.17) 


For small q , we recover S ~ q . Additionally, we can consider 
the case of finite q but small /?. For n 2 < 1, the RHS of 
Eq. (12.111) goes as 0 2 , whic h vanishes at 8 = 0 . This result 
was also found by iSrinivasan and Younj ( 20121 '). Therefore, 
at /3 = 0 any thin ring of an isotropic spectrum with k > q 
has no net effect on the zonal flow. On the other hand, for 
n 2 > 1 the 8 dependence in the RHS of Eq. ( 12 . 111 ) vanishes. 
Thus, at /3 = 0 a thin ring with k < q has a net damping 
effect on the zonal flow. 

For finite LS behaves a^f| 

S(x,n , m) = (4n 3 x) 1 [—n 2 (-l + m) + (1 + m) ^ — 1 — m 

+V / [(- 1 + n) 2 + m][( 1 + n ) 2 + m]j] + O (xf 3 ) • 

(2.18) 

For small /3, the /3 dependence cancels out of the RHS of 
Eq. CLED- Hence, zonostrophic instability is possible even 

2 Validity of this formula requires that 1 — n 2 is not too small. 

3 Validity requires that m ^ 0, because for m = 0 and n 2 < 1, 
the lowest order result vanishes. 


with P = 0. For concreteness, one might take m = 1, for 
which S simplifies to 

5(x,n, 1 ) = ^(- 1 + y7^)+0(x- 3 ). (2.19) 


Additionally, the small q limit of Eq. (I2T81) is 

n m 

S(x,n,m ) = 


X 2(1 + m) 2 


( 2 . 20 ) 


Thus, for an isotropic spectrum and finite LS goes as q 
at small q, rather than like g 4 as in the case of infinite L c [. 

These issues will be reexamined from another light in 
|Appcndix 2.4.A| where we give some physical understanding 
of the transfer of energy to long wavelengths. 

2.4-2.2 Instability of a Primary Wave to a Secondary 
Wave 

Zonostrophic instability can be understood in a very general 
way as the instability of some turbulent background spec¬ 
trum to a (zonally symmetric) coherent mode. As a special 
case, one can consider the background spectrum to consist of 
only a single mode. We show that in this case the dispersion 
relation of zonostrophic instability reduces exactly to that of 
the 4-mode modulational instability (sometimes called para- 
metric insta bi lity). Thi s corr espondence was first noted by 
ICarnevale and Martini dl982l ~) but was not understood in the 
context of the generation of zonal flows. 

The stability of a single, primary wave p to perturb ations 


1972; 

Gill, 

1974 

Gallagher et al. 


1974 1k rommesl. 120061 : IConnaughton et al 


201C 


2012). These calculations have used the 


fluctuating dynamical equations such as Eq. (ED and not 
a statistically averaged system. Generally one considers the 
unforced, undamped case, for which a single wave is an ex¬ 
act solution of the nonlinear dynamical equations. Concep¬ 
tually similar is the so-called secondary instability, where a 
growing, prim a ry eig e nmode gives r i se to a secondary mod e 


growing, primary eigenmode gives rise to a secondary mode 
( Rogers et all I2OO0I : IPlunkl . [20071 : IPueschel et all l2013l l. 


When the secondary mode grows much faster, the primary 
mode is treated as a stationary background. These sec¬ 
ondary instabilities can be more complicated, where, for ex¬ 
ample, the toroidal geometry of magnetically confined plas¬ 
mas results in the growing primary eigenmode having non¬ 
trivial spatial dependence. 

To calculate the stability of the primary wave using 
Eq. (12.11) . in general one needs to retain an infinite number 
of coupled, perturbing modes. However, typically one trun¬ 
cates the system, for example retaining a secondary mode q 
and the sideband pair p±q. Within this 4-mode approxima¬ 
tion and the further assumption that the primary has p y = 0 
such as a pure Rossby or drift wave and the secondary has 
q x = 0, the dispersio n relation for 4-mode mod ulational in¬ 
stability is given by ( Connaughton et al. I. l 20 ldl 


V 3 _ w 4 f 2M 2 (1 - s 2 )(l + s 2 + /)(! + /) 2 - (s 2 + /) 

V (l + /) 2 (l + s 2 + /) 2 ( S 2 + /) 

( 2 . 21 ) 

where X' = pX//3, s = q/p, f = p 2 L d 2 , M = ifiqp 3 f/3, and 
tfo is the amplitude of the background stream function. 
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Some studies investigated this phenomenon by using 
a form of CE2 where the inhomogeneity is assumed 
to vary slowly in spac e compared to the tu rb ulence 


(Dvachenko et all 1992; Manin and Nazarenkol, 

1994; 

Dubrulle and Nazarenko, 199?l; Smolvakov et al.. 

2000 b; 

Wordsworth, 

2009). With that assumption, the turbulence 


is described by a wave kinetic equation (see Section 5.1.1). 
The wave kinetic equation can also be recovered from the 
CE2 equation Eq. (I2.4al) by assuming dy -C d y , Taylor 
expanding, then Fourier transforming|j While those previ¬ 
ous studies are limited to the regime of small q, the CE2 
framework makes no assumption about the length scale 
of the inhomogeneity. Moreover, the previous studies did 
not draw a direct connection between the results from the 
statistical calculation and from the 4-mode calculation^ 

This dispersion relation Eq. ( 12.2111 can be recovered from 
CE2 and the zonostrophic instability dispersion relation 
Eq. (12.81) . To precisely compare, one must carefully select 
the background spectrum Wjj to correspond to a wave of 
stream function i/jq. If the initial background amplitude of 
mode p is ipo, then we write 

y) = </>o fe ipx ~ iut + (2.22) 

| Appendix 2.4.B| shows that this corresponds to a one-time, 
two-point covariance of streamfunction 

^ H {k x , k y ) = (27r)Vo [*(k - P) + <S(k + p)] . (2.23) 

From Eq. m , the corresponding covariance of vorticity is 
—4 

given by Wtf{k x , k y ) = k 'll n{k x , k y ), and thus, because of 
the delta functions, 

W H {k x , k v ) = (2-7t) 2 A [<5(k - p) + S (k + p)] (2.24) 

where we have defined A = ^{p 2 + L^ 2 ) 2 . There are two 
ways of achieving this background spectrum. First, we could 
choose the external forcing to be F( k) = 2 fiWfj. Since we 
want the dissipation term fi to disappear in the final expres¬ 
sion, fi can be chosen to be vanishingly small, in particular 
smaller than the eigenvalue A. Alternatively, as previously 
mentioned we could take the external forcing and the dissi¬ 
pation to be zero, in which case any arbitrary homogeneous 

4 The proper dependent variable to use for the wave action 
(see section 5.1.1.4.2) is A/"(k | y) = (1 — olzf k 2 L^ 2 )W (k | 

y ). For the CHME, this becomes J\T = ( k 2 /k 2 )W . With M as 
the dependent variable, the disparate-scale form of CE2 assumes 
wave-kinetic form. 

5 One reason a connecti on may not have be e n ma de is 
that the small-q result s in iManin and Nazarenko (1994) and 
ISmolvakov et alT d 2000 bl) based on the wave kinetic equation are 
incomplete. In their dissipationless (fi = 0) formulation, they ne¬ 
glect the term 2 if3qk x k y compared to A in the denominator of 
Eq. (12.91) . But this is invalid if A ~ q 2 because the neglected term 
is larger than the retained term. For example, when specialized 
to a single primary mode, both papers state that for the CHME, 
instability occurs when p 2 + L ^ 2 — 3p 2 > 0, and that \ ~ q 2 . One 
can obtain this result from the small q limit of Eq. (I2.25D if the 
/3 term is unjustifiably neglected. Careful analysis shows this re¬ 
sult also obtains, co rrectly, in the i/>q -» oo lim it. But contrary to 
statements made bv IConnaughton et al.1 d2010T) . the wave-kinetic 
formalism is not restricted to that large-amplitude regime. If the 
P term is retained, the full answer at small q can be recovered 
from the wave-kinetic formalism. 


spectrum trivially satisfies the CE2 equations. This latter 
point of view is closer to the traditional primary wave sta¬ 
bility calculations. 

Substituting Eq. (12.24|) into Eq. (|2.8|) . we find 


%X =2qApl 


P 


Py 


Py + 


hi 


Xpip 2 + 2 i/3qp x {p v - Xq) 


X P+P 2 + 2i/3qp x (py + X q ) 

(2.25) 


where dissipation has been neglected, pg_ = p 2 + ( p y ± q) 2 , 
and p± =p 2 ± + L~ 2 . 

When specialized to the case of a primary wave with p y = 
0 , the dispersion relation becomes 


- OnAr, 2 (l - t\ 1 2X P 2 +P 2 

q 2 9 PX \ p 2 ) 2 A2p4 p4 + 02g4p2 ' 


(2.26) 


Now, taking a zf — 1 and introducing the same normaliza¬ 
tions as used in Eq. (12.211) . we obtain 


s 2 + f 


A' = 


2As 2 X\l- s z )(l + s z + f) 


s‘ (0/p) 2 [ A' 2 (l + 6-2 + /) 2(1 + /) 2 + s 4] 

Letting A' = p 2 A//3 2 , after some simplification we find 
, / 4 ( 2^(1 — s 2 )(l + s 2 + /) — (s 2 + /) 


. (2.27) 


A ' 3 = AV 


(1 + /) 2(1 + S 2 + / ) 2 (s 2 + /) 


. (2.28) 


Since A' = p^ip q(1 + /) 2 //? 2 = M 2 (l + f) 2 , this exactly 
match es the dispersion relation given bv IConnaughton et, al.1 

(i 2010 ll in Eq. C2D- 

In the above calculation, we have shown that from CE2 
we recover the 4-wave modulational instability in the special 
case of a primary wave with p y = 0 and a secondary wave 
with q x = 0. In |Appcndix 2.4)Cl we generalize this to show 
that CE2 recovers the 4-wave modulational instability for an 
arbitrary primary wave and an arbitrary secondary wave. 

It may be at first surprising that the two dispersion rela¬ 
tions agree exactly, but retrospectively it makes sense. The 
4-wave modulational instability contains the primary wave 
p and the perturbations at wave vectors q and p ± q. From 
Eq. (12.8711 in |Appcndix 2.4.B| for the correlation between 
the primary mode k = px and sidebands k' = px ± qy, 
we see that the spatial dependence of the correlation goes 
as cos(pa; ± ^qy ± qy). Upon examining the CE2 calcula¬ 
tions, we see that the retained modes are the zonal flow 
8Ue± lqy (which corresponds to mode ±q) and the pertur¬ 
bations to the spectrum SW(k x , fcj / )e =l=z<?!/ . The perturbation 
5W (k x , k y ) is proportional to W}j{k x ,ky ± ^ q ), which is 
nonzero at k x = p and k y = i^q for the given primary 
mode. Therefore the perturbations kept within CE2 are pre¬ 
cisely the corresponding modes kept in the 4-mode trunca¬ 
tion. The CE2 instability calculation neglects higher har¬ 
monics of q such as e 2lqv at the linear level. These higher 
harmonics are precisely what is neglected by truncation to 
4 modes instead of retaining higher sidebands. 

In summary, the instability of a single primary mode can 
be thought of as a special case of the instability of an ar¬ 
bitrary background spectrum. In the fluctuating dynamical 
equations it is difficult to represent a turbulent spectrum as 
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an exact solution and hence calculate its stability. For this 
purpose, a statistical formulation such as CE2 is advanta¬ 
geous, since a homogeneous turbulent background can be 
represented more easily, particularly as a time-independent 
spectrum. 

When the homogeneous state is unstable, it is not obvious 
a priori what happens to the zonal flows. Do the zonal flows 
grow and saturate into a steady state? Or do they fluctu¬ 
ate turbulently, unable to persist in a steady state? Even 
though the answer is not obvious, within the Hasegawa- 
Mima or barotropic vorticity equation framework numerical 
simulations sometimes find steady zonal flows. But in other 
situations there may be fluctuating zonal flows. Complicated 
nonlinear physics determine what happens and it is difficult 
to determine what actually occurs without simulations. 

Another advantage of CE2 in particular is that it is possi¬ 
ble to calculate not only the instability of a turbulent back¬ 
ground, but also how the instability saturates. Analytic so¬ 
lutions are even possible in some regimes. This is undertaken 
in section l2.4.3.1l 


In addition, and especially relevant for this section, the 
CE2 simulations have yielded important information that in¬ 
form our analytic calculations. First, the simulations confirm 
that CE2, like both the quasilinear and original systems, 
exhibit zonal flows that can reach a steady state. Second, 
within the CE2 simulations one also sees the phenomenon 
of merging jets, which is u biquitous in DNS bu t has y et to 
be fully understood. Third. [Farrell and Ioannoul ( 20071 1 have 
found that CE2 can exhibit nonunique solutions, where the 
number of jets in the steady state depends on initial con¬ 
ditions. They also discovered that zonal flows emerge from 
homogeneous turbulence in a bifurcation triggered by zonos- 
trophic instability, and furthermore that the bifurcation is 
supercritical. Not only do these features guide our calcu¬ 
lations, they also demonstrate that CE2 possesses at least 
some of the essential physics of zonal flows as well as inter¬ 
esting and relevant nonlinear behavior. CE2 is therefore a 
system worthy of further understanding. 

In this second half of the section, we show that zonal flows 
can be understood as pattern formation. 


2.4.3 Pattern Formation 

In this first half of this article, we discussed the tendency for 
zonal flows to grow if they are not already present. Now, we 
consider what happens to such zonal flows beyond the initial 
stages of the instability. As the zonal flows grow larger, they 
reach an amplitude where they significantly modify the tur¬ 
bulence. Eventually, some kind of equilibrium between the 
turbulence and the zonal flows is reached. It is this saturated 
state that is of main interest in understanding the observable 
turbulence. 

Compared to zonal flow generation, the problem of zonal 
flow saturation has received much less attention in the the¬ 
oretical literature. Many of the works that have consid¬ 
ered it typically make an assumption of scale separation 
where the scale or wavelength of the zonal fl ows is much 
large r than the scale of the tu rbulence ( Diamond et all 


11998 : IConnaughton et all 12011 4. Zonal flows in plasmas, 
however, ar e often observed to be of comparable scale to the 
turbulence I Gupta et al.l . 120061 : iFuiisawa et al.l . [2004) . An¬ 


other line of inquiry is based on potential vorticity mixing 
(see Section 4.2). 

The CE2 equations provide a well-posed nonlinear system 
whereby the saturation of zonal flows can be investigated. 
As discussed at length previously, CE2 can describe the gen¬ 
eration of zonal flows through zonostrophic instability. But 
it can also describe the statistically steady inhomogeneous 
turbulence that results. 


(Farrell and Ioannou. 2003. 

20071. 2009|; Bakas and Ioannou 

2013bj: Oonstantinou et al. 

, 2013; Marston et all 2008 


2013) (see also 


sections 5.1.2, 5.2.2, and 5.2.4). Simulations of statistical 
equations, which evolve in time the covariances of the fluc¬ 
tuating fields, are distinct from conventional DNS, which 
evolve the amplitudes. The CE2 simulations have explored 
zonal flow physics in interesting ways distinct from DNS. 


2-4-3.1 Bifurcation of homogeneous turbulence 

Zonostrophic instability, which was discussed previously in 
Section 12.4.21 provides the starting point for our theoret¬ 
ical considerations. In some parameter regime, the coher¬ 
ent perturbations are stable and homogeneous turbulence 
is stable. But as a control parameter p, such as the fric¬ 
tion p or the forcing strength F, is adjusted, the homoge¬ 
neous state becomes unstab le ( Farrell and Ioannoul 120071 : 
ISrinivasan and Younel [20121 ). On either side of this instabil¬ 
ity threshold, or bifurcation point, the behavior of the sys¬ 
tem must be qualitatively different. Numerical simulations 
show that beyond the threshold the result is steady zonal 
flows and inhomogeneous turbulence. 

In mathematical terms, the statistical CE2 equations 
Eq. (12.41) possess translational symmetry y —> y + 5y. In 
other words, if {W(x,y \ y,t),U(y,t)} is a solution, then 
{W(x, y | y + Sy,t), U(y + Sy,t)} is too. Quite separately, 
when the system is zonostrophically stable, the homoge¬ 
neous solution, Eq. (12771) . manifests this symmetry by be¬ 
ing itself invariant to that transformation. When the con¬ 
trol parameter crosses the instability threshold, the system 
suddenly develops dependence on y as well as a mean held. 
The new solution is not invariant under translation. 

In order to fully understand the behavior of the system, 
analytic solutions would be beneficial in addition to nu¬ 
merical solutions. But the complexity of the nonlinear CE2 
equations means that finding a solution analytically is a 
formidable task and does not appear feasible in general. One 
way to proceed is by considering a regime where additional 
approximations can be made. Our approach is to investigate 
near the bifurcation point. The distance from the bifurca¬ 
tion point serves as a small parameter and facilitates further 
progress. 

The bifurcation analysis follows a standard procedure 
and involves a m ultiscale pertur bati on ex p ansio n about 
the threshold rtC ross and Hohenbergl . 1 19931 : iHovlel . 120061 : 
ICross and Greensidel . 20091 ). Since the bifurcation is super- 
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critical, only the lowest-order terms in the bifurcation anal¬ 
ysis are needed to provide saturation of the instability. The 
instability is known in the pattern formation literature as a 
Type I s instability. This type of bifurcation generically con¬ 
sists of a symmetry-breaking instability, a spatially periodic 
but temporally nonoscillatory marginal eigenvector, and a 
supercritical transition. 

We review the basic procedure of the perturbation expan¬ 
sion in an abstract n otati on. The full details are reported 
elsewhere ( Parker and Krommesi . l2013h b Consider a system 
with quadratic nonlinearity. Let (j> l> e an abstract vector, A 
be a linear operator, JV be a symmetric, bilinear operator, 
and F be external forcing. Let t = (p — pc)/pc be the nor¬ 
malized bifurcation parameter. Any of A, M, and F may 
depend explicitly on e. The basic equation is assumed to be 
given as 

0 = A </> + N(<j>, <)>) + F. (2.29) 

Given a nonzero equilibrium <j> e , we change variables by let¬ 
ting </> = <j> e + u to obtain 


0 = Lu + N(u, u ), 


(2.30) 


where Lu = Am + 2 N{<j> e ,u). In the context of the CE2 
calculation, (j> e = {Wjj, 0} and u = {W — Wjj, U}. By as¬ 
sumption, the linearization L around the equilibrium 4> e is 
stable for e < 0 , neutrally stable at e = 0 , and unstable for 
e > 0. 

The perturbation procedure employs slowly varying space 
and time scales in a multiple-scale expansion. We introduce 
the slow scales Y = e 1 / 2 ?/ and T = et, then let dy —> dy + 
and dt —> dt + tdy. Using these, we expand the 
operators L = Lq + C^ 2 L\ + eL 2 + c 3 ^ 2 Lq + • • • and N = 
Nq + e 1 / 2 ^! + • • •, and we expand the state vector u = 
e ' Mi + eM 2 + • • • . Collecting terms of the same order, we 
obtain the equations at O(e 1//2 ), O(e), and 0(e 3 / 2 ): 


< 1/2 ) : 

0 = Lqui, 

(2.31) 

Oil): 

0 = Lqu 2 + L 1 M 1 + Mo (mi, mi), 

(2.32) 

< 3/2 ) : 

0 = Lqus + L\u 2 + L 2 u\ 



+ 2M 0 (u\ , m 2 ) + Mi (mi, mi). 

(2.33) 


At O, Eq. (12.3111 states that mi is a null eigenvector 
of Lq, meaning it has a zero eigenvalue. The reality condition 
on m restricts the form of Mi to be 


mi = A(Y, T)r + A(Y,T)*r*, (2.34) 

where r ~ e lqcV and its complex conjugate r* are the right 
null eigenvectors, and A is the to-be-determined amplitude. 
These eigenvectors are periodic in y with critical wave num¬ 
ber q c , which is the first wave number to go unstable as e 
crosses zero. Once an inner product (•, •) is defined, then 
associated with the right null eigenvector r is a left null 
eigenvector l of Lq, which satisfies ( 1,Lqu ) = 0 for any u. 
The y dependence of l is also e lqcV . As is common in pertur¬ 
bative procedures, the amplitude A will be determined by 
nonlinearities occurring at higher order. 

At O(e), we first note that L\u\ = 0 automatically. This 
is because q c is marginally stable at the instability thresh¬ 
old: given a dispersion relation A (q, e) as a function of wave 


number q and control parameter e, then both \(q c ,0) = 0 
and d\/dq(q c ,0) = 0. The former equality yields Lqui = 0, 
while the latter equality yields the condition Limi = 0. In 
order to ensure that a solution for M2 exists, a solvability con¬ 
dition obtained by taking the inner product with the left null 
eigenvector must be satisfied. Applying this to Eq. (127321) . 
the solvability condition is 

(J,JVo(ui,mi)) =0 . (2.35) 

This solvability condition is automatically satisfied be¬ 

cause the quadratic nonlinearity implies Mo(mi,mi) ~ 1 or 
e ±2iq c y ^ w j 1 jj e 1 ^ g/lcy^ go the inner product (l, Nq(u\, mi)) 
vanishes. Therefore, given that a solution exists, one may 
write M2 as a linear combination of homogeneous and par¬ 
ticular solutions: 

U2=u 2 h+U2 P , (2.36) 

where 

u 2h = A 2 (Y, T)r + A 2 (Y , T)*r*, (2.37) 

Lou 2p = — Nq(ui, mi). (2.38) 

Since we have not yet determined A, we must proceed to 
higher order. Another unknown parameter A 2 has been in¬ 
troduced, but we will not need it in order to solve for A. 

At 0(e 3 / 2 ), note that L\u 2 h = 0 for the same reason 
that Liui = 0. Upon writing the solvability condition for 
Eq. (1031) . one finds that several terms vanish, leaving 

0 = (l, L 2 u\) + ( l , 21Vo(mi, u 2p )) ■ (2.39) 

This is the desired partial differential equation that deter¬ 
mines the amplitude A. It turns out that one never explicitly 
needs L\ or N\ in order to obtain this equation. 

We quote the results of the full analysis. After returning 
to the unsealed variables, the amplitude equation for A is 

c 0 d t A(y, t) = tc\A + C2 d^A - C3 1 A | 2 A, (2.40) 

where the Cj are order unity, real constants. All of the Cj 
should be positive (negative C3/C0 corresponds to subcritical 
rather than supercritical instability). 

Actually, one could have determined the form of 
Eq. (12.401) without going t hrough the actual calculation 
( Cross and Greensidel . I 2 OO 9 I '). The symmetries inherent in 
the original equation constrain the forms of possible terms. 
For instance, translational symmetry in y requires the am¬ 
plitude equation to be invariant to phase shifts of A, so that 
the lowest-order nonlinear term is uniquely determined to 
be \A\ 2 A. 

Furthermore, the behavior of Eq. (I2T01) is universal in 
the sense that, as long as all of the Ci > 0 , the qualitative 
behavior does not depend on the value of any of the Cj. 
This can be seen because a simple rescaling of A, y, and t 
eliminates Cj dependence from the equation. 

With Eq. (l2~40l) . the analogy between the zonal flows and 
the convection rolls in Rayleigh- Benard convection is com¬ 
plete. The transition to convection is governed by the same 
class of bifurcation and subject to the amplitude equation. 
The similarities between zonal flows and convection rolls al¬ 
luded to in section 12.4.1.11 are not merely descriptive, but 
mathematical as well. 
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The formulas for the Cj are complicated but are writ¬ 
ten in terms of the external parameters and integrals over 
the spectrum of the forcing. The formulas given in full by 
IParker and Krommesl 1 2013alb The important point is that 
it is possible to find a complete solution to the nonlinear 
CE2 equations, at least in a certain regime. 

There is an alternate method for obtaining Co, ci, and 
C 2 , which govern the linear behavior of Eq. (12.4011 for small 
A. The dispersion relation Eq. (12781) can be put into the 
form D( A, e,q) = 0 and can be Taylor expanded about the 
threshold. The conditions of the instability threshold require 
D( 0, 0, q c ) = 0 and dD/dq( 0, 0, q c ) = 0. Upon expanding D 
to lowest order about ( 0 , 0 , </c), one finds 

- 0, q c ) A = e^(0, 0, q c ) + 77 / 5 / 0 , 0, q c )(q - 9 c) 2 - 


dt 


2 dq 2 


(2.41) 

Up to a constant of proportionality, we identify cq = 
-dD/d\(0,0,q c ), ci = dD/de(0, 0, q c ), and C 2 = 
— ^d 2 D/dq 2 (0,0,q c ). This provides an independent check 
on the multiple-scale expansion calculation. However, this 
approach does not give C 3 ; for that one needs the full 
bifurcation calculation which includes nonlinear terms. 

Desired quantities of interest can be calculated analyti¬ 
cally from Eq. (12.401) . Linearizing about A = 0 and substi¬ 
tuting the form A ~ e* t e? ky , one calculates the growth rate 
to be 


A = 


eci — C2AU 

co 


(2.42) 


We recognize from Eq. (12.341) that k is the wave number rela¬ 
tive to q c so that k = q — q c - Steady state solutions including 
the nonlinear term also have the form A = A s (k)e lky , where 

,2 \ 1/2 


|A.(fc)| = 


eci — C2A7 
C3 


(2.43) 


The lowest-order correction to A is 0(e 2 ), while the lowest- 
order correction to A s is 0(e). 

Figure 12.41 verifies that Eq. (12.401) provides an adequate 
description of CE2 near the instability threshold. The ana¬ 
lytical growth rate found from Eq. (f2~42l) is compared with 
that from the exact dispersion relation Eq. Similarly, 

the analytical zonal flow amplitude found from Eq. (12.431) 
is compared with that from solving the full CE2 system as 
in section 12.4.41 We identify the amplitude A s of the first 
harmonic e lqcV with the numerically determined coefficient 
U\. The results are in excellent agreement. 

In addition to finding the steady states of the amplitude 
equation, one can ask whether those steady states are stable 
to small perturbations. Linear stability analysis about the 
solution A s e lky shows that it is unstab le to the E ckhau s 
instability when k 2 > eci/ 3 c 2 ( Cross and Greensidel . l2009l l. 

A stability diagram representing the various possibilities 
is shown in Figure 12.51 The neutral curve (N) indicates 
marginal stability of the A = 0 solution as a function of 
the wave number k and control parameter e. The A = 0 
solution is unstable to those k that are inside the neutral 
curve. At a fixed e > 0, steady-state solutions with A ^ 0 
exist at any of the k inside the neutral curve. The marginal 
stability of these A yf 0 solutions is indicated by the Eck- 
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Figure 2.4 Comparison showing agreement between numerical 
solution (blue circles) and analytic solution (black line), (a) 
Compensated growth rate A/e as a function of e at q = q c . (b) 
Growth rate A as a function of q at e = 0.01. (c) Compensated 
zonal flow amplitude t/i/e 1 / 2 as a function of e at q — q c . (d) 
Zonal flow amplitude Ui as a function of q at e = 0.0025. The 
zonal flow amplitude U\ is the first Fourier component of the 
zonal flow velocity U(y). (Adapted from New J. Phys. ©2014) 



Figure 2.5 Stability diagram for the amplitude equation 
EH). The labels ‘stable’ and ‘unstable’ refer to the nonzero-/! 
steady states. 


haus curve (E). Inside the E curve is a smaller band of wave 
numbers for which the steady-state solutions are stable. 

If a solution with an unstable wavelength is slightly per¬ 
turbed, it must evolve to reach a stable wavelength. The plot 
of Re A(y, t) in Figure 1231 with branches merging into wider 
branches, resembles similar plots of the zonal flow U ( y , t) 
in which jets merge. In the amplitude equation (12.4011 . the 
merging occurs in the nonlinear stage of the Eckhaus insta¬ 
bility. At the instant of merging ther e is a top ological defect 
known as a dislocation ( Cross and Greensidel . I 2 OO 9 I ). 

When the CE2 system is far from threshold, the amplitude 
equation ceases to be a quantitatively accurate description. 
However, many of the basic behaviors just described about 
the amplitude equation hold also for steady solutions of the 
CE2 equations, as we now verify by numerical solution. 
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Figure 2.6 Merging behavior in the amplitude equation it2~rni) 
[ReA(y, t) is shown]. (Prom New J. Phys. ©2014) 


2.4.4 Numerical Solution of CE2 and Stability 
Diagram 


In general, CE2 must be numerically solved. One approach 
is to evolve the CE2 equations in time until an equilibrium 
is reached. Our approach differs in that we solve the steady- 
state limit directly, i.e., Eq. (12.41) with d/dt = 0. We find 
steady-state solutions of zonal flows and turbulenc e using 
numerical techniques developed by Busse and Clever jBussel 
Il967l : Iciever and Buss j . IlDTdi : iBusse and Cleveil . Il 971)1) for 
the Rayleigh Benard convection problem. This method of 
solution uses a Galerkin expansion, where the dynamical 
variables are expanded in basis functions with unknown co¬ 
efficients and substituted into the equations of motion. The 
equations of motion are then projected onto the basis func¬ 
tions, yielding a set of nonlinear algebraic equations for the 
coefficients. 

The covariance of the turbulence and the zonal flow am¬ 
plitude are expanded as 

MNP 

W(x,y\y)= J2 J2 E Wmn P e imax e inb ye™, 
m=—M n=—N p——P 

(2.44) 

P 

U(y)= J2 u P eiP (2-45) 

p= — P 


where q is the fundamental wavenumber or 2n/q is the spa¬ 
tial periodicity of the zonal flows. There is a range of q that 
allows a solution. We obtain a system of nonlinear algebraic 
equations for the coefficients W mn p and Up by substituting 
the Galerkin series into Eq. (12.41) and projecting onto the ba¬ 
sis functions. To demonstrate the projection for Eq. (l2~4al) . 
let <t>mnp = e lrnax e mbv e xpqv . We project Eq. (!2.4aP onto 4> rs t 


by operating with 


/ 27T 27T 27t\ 
\ a b q ) 


— 1 rTc/a n'K/b f‘' 7T /Q 

/ dx dy dy first- 

J —7T /a J—n/b J—Tr/q 


r v/b 
f —n/b 


n/q 

-n/q 


(2.46) 


Projection of the first term, ( U+ — U—)d x W , yields 
Irstp'mnpUp’Wmnp , where repeated indices are summed 
Over, Irstp'mnp itYKlSm.. r ftp' -j-p—t.O (uy rr_), (Ty — 

sinc(a± 7 r/b), and ce± = nb — sb ± \p> q. The other terms 
of Eq. (12.4al) . as well as Eq. (I2.4bl) . are handled similarly. 
In total, we generate as many equations as there are coeffi¬ 
cients. 


The system of nonlinear algebraic equations is solved with 
Newton’s method. One feature of Newton’s method is that 
it requires a good initial guess. We attain a suitable guess 
by using the bifurcation solution near threshold. Then we 
adjust a parameter in small increments towards the desired 
value, a technique known as numerical continuation. The 
solution at the previous value of the parameter can serve as 
the initial guess for the next value. 

Once a steady-state solution is found, its stability can be 
assessed. The general method, again following Busse and 
Clever, involves linearizing the equations of motion about 
the steady state. Since the equilibrium is periodic in y, the 
perturbation may be expressed as a Bloch state. Then the 
perturbation is expanded in the same Fourier-Galerkin basis 
functions used to express the equilibrium. The equilibrium is 
unsta ble if there are any eig envalues with positive real part 
( Cross and Greensidel . 120091 ). Further details of the numer- 
ics rega rding the equi l ibrium and stability may be found in 
IParker and Krommesl 1 2013bl ial). 

In the same manner as for the amplitude equation, the 
results are organized into a stability diagram. Figure 12.71 
displays the stability diagram for the CE2 system with 
infinite deformation radius. The control parameter on 
the y axis is 7 = e 1 / 4 /! 1 / 2 /r -5 / 4 , a fundamental param¬ 
eter c ontrolling the je t dyn a mics ([ Danilov and Guraric, 


2004 ; 


Gnlj 3 erm = etal 1 1201(1 IScott and Dritschel 


2012; 


Tobias and Marstonl . 120131 : iBouchet et al. . 20131 ). This pa¬ 


rameter is related to the zonostrophy parameter Rp by 
7 = R^j. Near the instability threshold, the stability dia¬ 
gram resembles that for the amplitude equation (see Figure 
1231) . as it should. The Eckhaus (E) instability forms the 
stability boundaries near the threshold in the sense that if 
one starts inside the stable region and increases or decreases 
q, the Eckhaus instability is the first instability triggered 
unstable. Farther from threshold, at larger 7 , other instabil¬ 
ities form the boundary (Li and Ri in the diagram). These 
other instabilities have not yet been studied in detail. 

In Figure 1.2 of Section 5.2.2, Farrell and Ioannou show 
a similar stability diagram for the j3 plane. Their statisti¬ 
cal approach, called S3T, is mathematically equivalent to 
CE2 although a different coordinate system and numerical 
method are used in practice in the computations. Unlike our 
numerical method, which develops problems at larger values 
of 7 , their method has no problem achieving values of 7 far 
from the critical value. In their figure, as the strength of 
the forcing is increased (which corresponds to increasing 7 ) 
well beyond the critical value, the region of stability curves 
to the left toward small wavenumbers or larger jets. Such 
is the behavior qualitatively expected in order to follow the 
Rhines scaling. 


2.4.5 Summary 

I 11 the first part of this article, we joined numerous other 
authors in offering a perspective on the generation of zonal 
flows. We found a deep connection between the stability of a 
single wave and the zonostrophic instability of homogeneous 
turbulence. In particular, the 4-wave modulational instabil¬ 
ity can be recovered exactly as a special case of zonostrophic 



















































14 



Wave Number q 


Figure 2.7 Stability diagram for the CE2 equations. For 7 
above the bottom of the neutral curve (N), the homogeneous 
turbulent state is zonostrophically unstable and the result is 
inhomogeneous turbulence with zonal flows . Ideal states are 
stable within the marginal stability curves E, Li, Ri. The 
stability curve is consistent with the dominant zonal flow wave 
number from independent QL simulations (crosses). The 
stationary ideal states vanish to the left of curve D. The black 
dashed line depicts the Rhines wave number. (Adapted from 
New J. Phys. ©2014) 


instability within the CE2 formalism. In addition to a sin¬ 
gle wave, we also examined the case when the background 
spectrum is isotropic. When the deformation radius is fi¬ 
nite, there are some notable differences in the physics of 
eddy forcing of zonal flows, especially for long-wavelength 
jets. 

In the second part, we considered zonal flows beyond the 
initial stages of growth and asked how they saturated into 
a steady state. We described zonal flows as pattern forma¬ 
tion amid a bath of turbulence. A deep understanding of 
the spontaneous symmetry breaking of statistical homogene¬ 
ity attained through the CE2 framework reveals behaviors 
such as the existence of multiple solutions with different jet 
wavelengths and the phenomenon of jet merging to reach 
a stable wavelength. These features have been observed in 
simulations. 

The pattern formation view of zonal flows is quite gen¬ 
eral. It possesses a far broader scope than the minimalistic 
2D models considered here. The behaviors predicted by the 
amplitude equation should be expected any time there is 
a spontaneous symmetry breaking with the appearance of 
steady zonal flows. For example, in a generalization of the 
Hasegawa-Mima equation that includes a resistive instabil¬ 
ity, so me o f th e exp ected features occur along with zonal 
flows ( Numata et all 12007 ). 

We have been emphasizing the role of symmetry break¬ 
ing. But in reality, a 0 plane does not exist. Moving to a 
more physical model such as the surface of a rotating sphere 


destroys the north-south translational symmetries associ¬ 
ated with a ft plane. Do any of these results apply to zonal 
flows in spherical geometry? Although this question should 
be studied in detail, we offer one possibility. Due to the 
latitudinal variation of the Coriolis parameter, the turbu¬ 
lence is always inhomogeneous on the sphere. A transition 
from homogeneous to inhomogeneous turbulence is not the 
right description, but perhaps some type of transition may 
still occur. Besides for the development of inhomogeneity, 
another aspect of the bifurcation on a ft plane is the spon¬ 
taneous formation of a mean held, i.e., the zonal flow. We 
suggest that this mean-held generation may survive for flow 
on a rotating sphere, and would be observable as a control 
parameter is varied. The zonal flow still behaves as an order 
parameter in this more general type of scenario. 


Appendix 2.4.A Zonostrophic instability and the 
physics of disparate-scale interactions 


We found in Sec. 12.4.2.11 using the CE2 approximation, that 
for the zonostrophic instability the behavior of the effective 
forcing on the zonal flows depended on whether the defor¬ 
mation radius L^ was infinite or finite. For L^ = 00 , wave 
numbers k > q of an isotropic spectrum produce no net forc¬ 
ing, whereas there is net forcing for finite L^. We also saw 
in Sec. 12.4.2.21 that the standard equations for generalized 
modulational instability are a special case of those for zonos¬ 
trophic instability. In this appendix we will discuss some of 
the connections between these various results. 

Various proposed mechanis ms for the formation of zona l 
jets have been summarized bv lBakas and Ioannoul ( 2013al ). 
They listed “turbulent cascades, modulational instability, 
mixing of potential vorticity, and statistical theories”; their 
work focused on the implications of the S3T closure. As they 
pointed out, one of the key points to be reckoned with is 
that “previous studies have shown that shearing of isotropic 
eddies on an infinite domain and in the absence of dissi- 
pation and 0 does not produce any net mom entum fluxes 
( Shepherd! . 1 19851 : iFarrell [~1987l : iHollowavL I 2 OI 0 I ) .” Note that 
none of Bakas and Ioannou, Shepherd, or Farrell cited the 
closely rela ted, detailed, and compelling discussion given by 
iKraichnanl ( 19761 . Sec. 5) of the physical mechanisms that 
underlie long-wavelength flow generation for 2D Navier- 
Stokes turbulence in both coherent and stochastic situa¬ 
tions. The implications of that work also do not seem to 
be ap preciated by many workers on the modulational insta- 


be ap preciated by many workers on the modulational insta¬ 
bility. |Holloway| ( 201(3) did cite it, discussed why the works of 


Kraichnan and Shepherd seem to have had limited impact, 
and went on to provide valuable new insights about some 
of the apparent contradictions that arise in various descrip¬ 
tions of eddy shearing. Our discussion below adds additional 
perspectives. 

Although Holloway provided some description of Kraich- 
nan’s calculations, we find it necessary to discuss them here 
as well. (Essential background can be found in the article 
by Krommes and Parker in Sec. 11.11 of this book, which we 
will abbreviate as KP.) Kraichnan’s original analysis was 
for 2D homogeneous Navier-Stokes turbulence (for which 
Ed is infinite). The analysis, which is generalized here to fi- 
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nite Ld, turns out to be useful not only for understanding 
long-wavelength flow generation in homogeneous turbulence 
(see KP), but also for gaining an intuitive understanding 
of the physics of zonostrophic instability and bifurcation to 
inhomogeneous turbulence. We find a connection to various 
limits reported in Sec. 12.4.2.11 In order to provide necessary 
background, we will first review Kraichnan’s original anal¬ 
ysis; we will also discuss how it is related to conventional 
calculations of modulational instability, thereby making a 
connection to our observation in Sec. 12.4.2.21 that modula¬ 
tional instability is a special case of the zonostrophic insta¬ 
bility. Then we will generalize the basic ideas to situations 
with finite L^. For those cases, we will show that some of 
Kraichnan’s conclusions are nontrivially modified in a way 
that is consistent with the results described in Sec. 12 . 4.2.11 
and we will provide some heuristic understanding. 


2-4-A.l Review of Kraichnan’s discussion of negative 
eddy viscosity 


Kraichnan framed his analysis as a calculation of a sta¬ 
tistical eddy viscosity p(q | fe m ; n ) felt by resolved scales 
(wave number < k m j n ) due to the interactions with un¬ 
resolved sub-grid scales; see the discussion in Sec. 5.1.4.2 
of KP. In the asymptotic limit q -C fc m ; n , he found that 
p(q I L m) < 0 in 2D; this is th e famous negative eddy 
viscosity. Krommes and Kiml ( 2000l ) discussed an important 
connection between that result and a certain formula for 
the rate of zonal flow generation, and aspects of that anal¬ 
ysis will be useful here as well. Kraichnan also pointed out 
that p(q | fc m in) actually vanishes for situations in which 
the interactions are dominated by long-wavelength strain¬ 
ing of turbulent excitations confined to k > fe m i n - We will 
generalize that result to models with finite L^. 

Kraichnan described the transfer of energy from short to 
long wavelengths in 2D turbulence as resulting from the gen¬ 
eration of a ‘secondary flow,’ a concept closely related to the 
mechanis m of ‘secondary ins t ability’ consi d ered by various 
authors ( Rogers et al ], [2000 ; El un 2, 120071 : Pueschel et all 
1201 dh . He began with a blob of short-wavelength vorticity 
(having central wave vector K = Ky) initially localized 
within a circular domain of radius D (KD 1) and pos¬ 
sessing the stream functior@ 




( Kquq 


f(x) cos (K ■ x), 


(2.47) 


where A'o = K (0) and 


f(x) = exp 



where p 2 = x 2 + y 2 . (2.48) 


6 Kraichnan used Kq = 1 and uq = 1, but we prefer to leave 
them general so that the dimensions of various quantities are cor¬ 
rect. We have changed some of his notation as well. For example, 
we have used uppercase K and Q for the specific wave vectors of 
the turbulence and the straining field, respectively. 



-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 


Figure 2.8 The velocity field corresponding to Eq. 112.471) . 
Lengths are normalized to D\ K y = 27r. 



-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 

Figure 2.9 The vorticity field corresponding to Eq. i2~m . 


The resulting velocity is 

u = z X V?/> (2.49a) 

= - (Jff'j “0 f(x)[z X K sin (IT • x) 

+ (KD)~ 1 z x (x/D) cos {K ■ *)] (2.49b) 

(Fig. 12.81) . and the vorticity is 
u: = V 2 -!/’ (2.50a) 

= -K 0 u 0 f(x)({[ 1 + (KD )~ 2 {2 - (p/D) 2 ]}cos(K ■ x) 
-2(KD)~ 1 K- (x/D)sin(K -x)) (2.50b) 

(Fig. 12.911 . 

One way of understanding the role of the shaping func¬ 
tion / is by inquiring about the spectral content of ip- One 
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has 



f(x) cos (K ■ x ) 

2 + e - i \ k + K \ 2D2 


(2.51a) 


(2.51b) 


With Q = D _1 (underlying vector Q 's will be intro¬ 
duced later), this describes a spectrum containing the pri¬ 
mary mode K and all sidebands having magnitudes up to 
p± = \ t k-q\. The role of / is thus to introduce side¬ 
bands that are necessary in order that triad interactions 
can occur between the primary, the sidebands, and a long- 
wavelength disturbance with characteristic wave vector Q ; 
compare the minimal system of four wave vectors K, P±, 
and Q [KP, Fig. 5.4 (right)] used in modulational-instability 
calculations. 

Kraichnan now introduces a long-wavelength straining 
field having potential 


///// / HU \ \\\\ 

yy///n n \\\\\y 

//////{ i \\\\\\\ 
x / / / t \ \ \ \ \ v' 




1 V \ \ \ x 

t v \ v *. --k. 


\ \ i / / ** . 

x, \ \ \ t t / / s s’. 

X \ \ \ \ t t t / /S. 

\\\\X\\\!i f //// 
V\W\Ht t t f //XX 
\\\\\Ht ft ///// 



Figure 2.10 The long-wavelength straining field. 


V(x) = —axy, (2.52) 

where a is an unspecified amplitude (having the dimensions 
of frequency or vorticity). The straining velocity is 

v(x) = zx VV r = a(x x — yy) (2.53) 

and is visualized in Fig. 12.101 This is a flow with pure rate 
of strain, i.e., it is irrotational: z ■ V X v = V 2 F = 0. To 
make contact with the calculations of modulational instabil¬ 
ity, consider its spectral content, which is 

V(q) = (2n) 2 a$'(q x )$ l (q y ) = (2tr) 2 aS' (q). (2.54) 

This is a somewhat unusual and pathological function. How¬ 
ever, it can be regularized by replacing the derivatives of the 
delta functions with finite-difference representations, e.g., 
5 1 (q) ~ [ S(q + Q) — 5(q — Q)]/2Q for small Q. One is led 
naturally to this approximation by noting that since v(x ) 
will be interacting with the shaped blob of short-wavelength 
vorticity, which localizes distances to Qx < 1, it makes little 
qualitative difference if one replaces V(x) by 

V(x) = —aD 2 sm(Qx)sm(Qy) = V+(x) — V-{x ), (2.55) 

where 

V± = ^aP 2 cos(Q± ■ x) (2.56) 

with 

Q ± =Q(S±y). (2.57) 

One has 

V±(q) = n 2 aD 2 [S(q - Q±) + 5(q + Q±)\. (2.58) 

Thus the original irrotational straining field is the differ¬ 
ence of two fields, each possessing both strain and vortic- 



Figure 2.11 The velocity field corresponding to V+. It is built 
from Q .|_ = (1, 1) T and has both strain and vorticity. 

ityd whose wave vectors Q± are oriented along the ±45° 
diagonals, as illustrated in Figs. [2J~T1 and [ 2 ~. 121 

Conventional modulati onal instability calculations 
( Nazarenko et all . ISec. 4,d . and references therein) begin 
with a single Q and its negative. The initial state of the 
instability thus possesses both vorticity and strain. We will 

7 In plasma physics and possibly elsewhere, it is ubiquitous to 
illustrate physics related to eddy ‘shearing’ with velocity fields 
like the one shown in Fig. 12.111 which possess vorticity as well 
as strain. Usually the rotational part of the interaction is not re¬ 
marked upon. While that often does not matter for simple heuris¬ 
tics, some arguments and illustrations would be clearer if Kraich- 
nan’s example were followed and a field with pure rate of strain 
were used. 
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Figure 2.12 The velocity field corresponding to V —, built from 
Q~ =(1,-1) T . 


Figure 2.13 The secondary flow, containing four vortices, that 
arises by self-interaction of the small-scale motion. 


see shortly how such an instability is related to Kraichnan’s 
procedure. 

We continue to review his calculations. The next step 
is to find an expression for the time rate of change of 
short-scale energy due to the straining. The vorticity equa¬ 
tion dtu + v ■ Vo; becomes, in a Lagrangian representa¬ 
tion, cLo/dt = 0 with the characteristic equations dx/dt = 
v(x) = a(x, —y) T - As Kraichnan observed, it follows that 
an initially circular blob is stretched into an ellipse with 
major axis in the x direction, while the central wave vec¬ 
tor is stretched according to K x {t ) = K x ( 0)e -ot , K y (t) = 
Ky{0)e at . Direct calculation of the time rate of change of the 
spatially-integrated energj@ £ = ^u 2 oc K~ 2 + 0((KD)~ 2 ) 
then leads, for K( 0) = K y y, to the initial energy loss rate 
£ = —2 a£ to lowest order. By considering the secondary flow 
that is generated by u (i.e., by evaluating dtAuj = —u ■ Voj) 
at t = 0), Kraichnan demonstrated that the lost energy 
shows up as energy of interaction between the secondary flow 
and the straining flow. Figure d. 13l illustrates that secondary 
flow, which is such as to reinforce the original straining flow 
near the origin (for positive a). 

This nonrandom mechanism, with energy transfer medi¬ 
ated by the amplitude a, is closely related to conventional 
calculations of modulational instability. Those describe an 
eigenvalue problem in which the unstable eigenvector pos¬ 
sesses both strain and vorticity and grows coherently. In 
Kraichnan’s calculation, the original straining field is rein¬ 
forced by the vorticity of the secondary flow. If that rein¬ 
forced field were taken as a new initial condition and the pro¬ 
cess were repeated, the evolving long-wavelength flow would 
contain vorticity as well as strain, as in the modulational- 
instability calculations. To understand the time scale for the 


reinforcement, consider the secondary-flow equation 

dtAuj = —u ■ Vw. (2.59) 

It is straightforward to use the results (I2.49bl) and (I2.50bl) 
to show that the secular part of the right-hand side of 
Eq. (12.5911 is at t = 0 

(U ■ 

= -2uo/ 2 D _4 (z • Kx x)(K-x) (2.60a) 

oc (g) (J) fiQuof ~ (Quo) 2 . (2.60b) 

The frequency Quq is the circulation rate or vorticity of 
one of the vortices shown in Fig. 12.131 That should also 
be the characteristic rate of the reinforcement, so one con¬ 
cludes that the characteristic rate for the growth of the long- 
wavelength flow is A ~ Quo . This agrees with the result of 
a modulational-instability calculation in which wave effects 
are neglected and the asymptotic limit of small Q/K is taken 
( Krommej . boQfil 'l; it is also consistent with the implications 
of Eq. KZ?m . 

We now turn to the implications of this analysis of co¬ 
herent interactions for statistical scenarios. Kraichnan ad¬ 
dressed this in the context of the 2D Navier-Stokes equa¬ 
tion; we are interested in the generalization of his analysis 
for cases with finite deformation radius. The basic calcula¬ 
tion makes the straining amplitude a a random function a(f) 
and also assumes that the wave vector K of the small-scale 
motion is oriented randomly, having angle <j> with respect to 
the y axis. With the assumption of passive advection of the 
small scales by the straining field, it is straightforward to 
find that 

K 2 (t)/K 2 = cosh[26(t)] + sinh[26(t)] cos(2c)>), (2.61) 


where 

8 There is a crucial misprint in Kraichnan’s formula for the 
initial kinetic energy in the second line after his Eq. (5.7); a factor 
of k ~ is omitted. 


b(t) = f dta(t). 

■ ' 0 


(2.62) 
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Kraichnan noted that (K 2 )^ (averaged over <j> but not a) 
typicall y grows in mean square, consistent with general re¬ 
sults of Cockel ( 19691 '). However, according to Eq. (I2.49bl) the 
2D Navier-Stokes energy £ = ^u 2 is proportional to A' -2 , 
and Kraichnan found that 


<A- 2 (t )) 0 = A 0 - 2 , (2.63) 

i.e., £ is independent of time in spite of the random stretch¬ 
ing and squeezing of K(t). It is worth quoting Kraichnan’s 
interpretation of this in his own words, since we will shortly 
give a more general discussion. He was concerned with the 
physics of the isotropic 2D eddy viscosity, which we repeat 
here for conveniencaj: 


I K 


f 

Jk„ 


dkd, 


qkk 


d[k 2 U(k)} 


dk 


(2.64) 


Regarding Eq. (12.6411 . he observecF^l 


reader would be forgiven for pondering whether it is energy 
diffusion that is meant, but further thought and reference 
to the discussion of nonlinear invariants in KP, Sec. 1.1.4.1, 
lead one to conclude that it is actually enstrophy that dif¬ 
fuses (total enstrophy being the nonlinear conserved quan¬ 
tity). We will see that a proper understanding of this latter 
point will also clarify the role of 9 qkk -, it is the autocorre¬ 
lation time associated with the wave-number diffusion co¬ 
efficient D k of enstrophy, which is more fundamental than 
6 itself. 


2.4-A.2 The effects of finite deformation radius 

It is instructive to consider these issues for cases involving 
finite deformation radius L specifically the Charney- 
Hasegawa-Mima equation (CHME) and the modified 
Hasegawa-Mima equation (mHME). 


“The integrand is a total derivative except for the k dependence 
of 9kkq • This means that any addition to the spectrum U(k) for 
k > A: ln j n which vanishes at k = A: lmtl would add nothing to 
ji(q | fc min ) were it not for the k dependence of 9 kkq .” 

He then interpreted the results of his model calculation as 
follows: 

“If 9kkq is dominated by low-wavenumber straining, in correspon¬ 
dence to our present discussion, it is independent of k and the 
integrand of [Eq. (12. 64 11 1 is a total derivative. Thus any excita¬ 
tion, described by Vt(k), which is totally confined tcFl k > /C m in, 
gives zero contribution to the effective eddy viscosity exerted on 
q <C k. This is a direct consequence of [Eq. I12.63D 1 which says 
that low-wavenumber straining of the small scales gives a diffu¬ 
sion process in wavenumber with no average loss of kinetic energy. 
By conservation, there is then no net gain of kinetic energy by 
the straining scales. On the other hand, if k m i n falls within the 
small-scale excitation, the diffusion of the excitation to smaller k 
occurs at wavenumbers < k m [ n and is not counted in [Eq. (12.64D 1 
which then includes only the outward diffusion. The latter does 
involve a net loss of kinetic energy by the small scales and thus 
gives rise to a negative contribution to the eddy viscosity.” 

Kraichnan’s insights here are deep and important, but 
two points require further discussion. First, he attributes 
the nonvanishing of /i q to the k dependence of Q q kk-> but he 
does not give a satisfactory explanation of why that quantity 
should be fundamental. Second, he uses the phrase “diffusion 
process in wavenumber” without clearly specifying exactly 
what quantity is diffusing. Given the immediate context, the 

9 Following the conventions used by KP, we indicate discrete 
Fourier transforms by subscripts (e.g., A/fc) and integral trans¬ 
forms by arguments [e.g., A/”(fe)]. Two-point spectra are nor¬ 
malized such that the fluctuation intensity is J\f = yfy A/fe = 
(27r) -d jdkj\f{k ), where d is the dimensionality of space (= 2 for 
the present discussion). The velocity spectrum is tiy.. Energy and 
enstrophy spectra are defined with a factor of ^ relative to IA. 
We will not have occasion to use omnidirectional spectra such as 
the common E(k ), which incorporate the wave-number volume 
element. 

10 For consistency with our notation and that of Kraichnan’s 
model, we have interchanged k and q from Kraichnan’s original 
usage in his Sec. 4. We also write k m instead of km and U instead 
of U. 

11 The published text contains the typographical error k < 
fc m in instead of the correct k > k m i n . 


General formulas for energy gain and loss — The rel¬ 
evant nonlinear invariant is (see the background material in 
KP, Sec. 1.1.4.2) A/fc = where cr\ = k 2 for the CHME 

and a k = k for the mHME. Here k = a k + k , where 
a k = 0 for zonal modes and a k = fc 2 otherwise ( k^ = Lf 1 )-, 

also, £ k = \k 2 (\S(p k \ 2 ). (The 2D Navier-Stokes case is re¬ 
covered for a k = 0.) We assume a h omogeneous ensem¬ 
ble wi th random long-wavelength flows. iKrommes and Kiml 
f : 200 oh showed that, upon expansion in t = q/k -C 1 of an 
anisotropic extension of Kraichnan’s test field model, a dif¬ 
fusion equation ensues for the short-wavelength spectrum: 


dt dk k dk 


(2.65) 


where 


D fe = 2A: 2 (4) ^2(QQ)\kXq\ 2 (j^J d k ,_ k , q Af<. 

( 2 . 66 ) 

(Krommes and Kim also gave a heuristic random-walk 
derivation of D^,.) Short-wavelength energy £ k evolves ac¬ 
cording to the nonconservative equation 


9i = ii.n d( g fc g fc) 

dt <j| dk k dk 


(2.67) 


Upon writing this as much as possible in conservative form, 
one finds 


dt 


d_ 

dk 


Dfc 


ML 

dk 


JL 

dk 


-2 


2D fc 


dk 


+ dk' Dk ' 




dk 


a/T 


( 2 . 68 ) 


Thus, while wave-number diffusion (first term) does act on 
the short-scale energy, £ k also experiences drag (second 
term) as well as an intrinsic loss mechanism (last term). 
The loss term describes the second-order, statistically aver¬ 
aged effect of random refraction of the ray trajectories of 
the small-scale wave packets; it is built from the first-order 
refraction effect discussed by KP, Eqs. (5.82) and (5.83). 
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Mathematically, it arises because the scale factor a k that re¬ 
lates A/fc and £ k does not commute with the Poisson bracket, 
involving large scale X and large wave number fc, that gen¬ 
erates weakly inhomogeneous wave kinetics. 

To verify that energy lost from the short scales shows up 
in the large scales, consider the e quation for long-wavelength 


m tne large scales, consider tnc equation tor iorig-wavcicngtn 
energy evolution £ q , which from iKrommes and Kiml ( 200nh 


9f£q — 2 'JqSq 


with 

7 q 


(2.69) 

0M> 


= ( = 2 ) E £ (=2 ) |feXq | 2 k qe q ^ k , k q 


dk 
(2.70) 

It is then straightforward to verify the energy conservation 
law 


d t £' 


= E ^( 27,0 = -a t e > = -^-LJL.D k . ( - 

u k 


9^1 

dk 


(2.71) 

by integrating the last expression by parts. We ignore sur¬ 
face terms, meaning that we consider excitations entirely 
localized within the domain of integration. 

The form of Eq. (12.681) can be used to give further insight 
to Kraichnan’s observation that in the isotropic 2D Navier- 
Stokes case the energy transfer would vanish for localized 
excitations were it not for the k dependence of 9 qkk . Clearly 
the first two terms contribute nothing; they merely rear¬ 
range short-scale energy locally in k space. The last term of 
Eq. (12.681) can be written as 


{d k 


r kWk 


(2.72) 


where the ‘flux of inverse scale factor’ is 


Tfe = ~D fe • dkcr k 


(2.73) 


It is this term, the statistical manifestation of the ray equa¬ 
tion k = — VOfc, where Q k is the nonlinear advection fre¬ 
quency (see the discussion of the first-order distension rate 
yjJ by KP, Sec. 5.1.4.2), that has the potential to trans¬ 
fer energy. Because of the factor of J\f k in Eq. (12.721) . the 
term is not conservative. Rather than describing a rate of 
redistribution of energy among the small scales, T*. gives 
the rate of transfer to the secondary flow and thus to the 
large scales. But if the divergence of that flux vanishes, 
no net energy transfer ensues (no secular contributions to 
secondary flow are generated). There are two contributions 
to that divergence, namely the k dependencies of D*. [oc 
k2 (°-l/k 2 ) 2 9 q ,k-k] andof d fc o-~ 2 = -cr“ 4 S fc cr| = -2<r“ 4 fc. 
Note that the dependence on a k cancels out between this 
term and D^. For the isotropic 2D Navier-Stokes case 
(a k = 0), one finds T fc oc k(k 2 9 qkk ) x ( k ~ 3 ) = kk~ x 9 qkk \ 
thus d k • Tfc vanishes to the extent that 9 qkk is independent 
of k. For the cases with finite deformation radius, the re- 

^' 1 A ^ 

suit is instead T*. oc kk~ (fc /fc )9 qkk , which has nontrivial 
divergence even if 9 qkk is independent of k. One sees that 
Kraichnan’s result that the energy transfer is controlled by 
the k dependence of 9 qkk is a special case; of more funda¬ 


mental relevance is the k dependence of which stems 
from the underlying physics of random ray refraction. 

Upon summing Eq. (12.681) over fc, one finds that the ex¬ 
plicit result for a localized isotropic spectrum is 



(2.74) 


By virtue of energy conservation, this reduces to 
— 2 ^ a 7 qfq where 7 q = —q 2 fi q and y q , which general¬ 
izes the 2D Navier-Stokes result (12.641) . is 


I fc n 


n OO 

" &mir 


dk 9k Kfv H 

(2.75) 

For otk = 0 (q = q and k = fc), this reduces correctly to 
Kraichnan’s result (12.641) . 

The interpretation of the ratio Rk = fc 2 /o | is that it is a 
measure of the portion of the physics devoted to perpendic¬ 
ular advection. To be specific, we discuss the plasma case. 
The Hasegawa-Mima equation for the magnetized plasma, 
Eq. (V.1.49), embodies the two quite different physical pro¬ 
cesses of (i) perpendicular advection of vorticity (the Vjj 
term), and (ii) parallel electron response, which is rapid and 
adjusts essentially instantaneously to changes in the electro¬ 
static potential (the a*, term). The total energy is the sum 
of (i) the kinetic energy associated with the perpendicular 
flow, and (ii) the compressional energy associated with the 
parallel motion. R k is the fraction of total energy associ¬ 
ated with the perpendicular processes. (It approaches 1 for 
a mode whose wavelength is much smaller than L ^.) It is 
only that fraction that is relevant for the random ray refrac¬ 
tion. More directly, the presence of Rk in Eq. (12.661) for Dj, 
arises from the fact that the effective frequency for advec¬ 
tion of the short scales is reduced for the CHME by a factor 
of Rk from the nominal k ■ V q of the mHME; it appears 
squared because the random nature of the refraction leads 
to wave-number diffusion, a second-order effect. 


Generalization of Kraichnan’s model to include fi¬ 
nite deformation radius — We now show that these 
results are consistent with a generalization of Kraichnan’s 
model. For definiteness, we consider the modified Hasegawa- 
Mima equation. I 11 order to construct a stream function that 
corresponds to a short-scale blob of generalized vorticity, 
and in view of the forms of £k and A/fc, one must replace A' 2 
in Eq. (12.471) by A' 2 (but not K by K). Because £k = k “A fk 
and N is conserved under the disparate-scale interaction, it 
is useful to consider 

R(t) = I<i(K~ 2 (t )>* (2.76a) 

1 r 2 ' K k 2 

= E/ d<t> -=-- 

2n Jo ax + [cosh(26) + sinh(26) cos(20)] A'g 

(2.76b) 

= {1 + 2ocosh[2&(t)] + a 2 }” 1 / 2 , (2.76c) 
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where a = a/K q. (ajc = k^ 2 for short-scale modes.) This 
correctly reduces to Kraichnan’s result (12.6311 for a = 0, but 
depends on the random straining otherwise. 

At t = 0, one finds 


R( o) 


An 


a + A'q 



(2.77) 


This is trivial (straining has not yet acted at t = 0); it should 
not be confused with Eq. (I2fo3ll . which holds for all times, 
and merely confirms that an average over an isotropic wave- 
number distribution does not change the isotropic quan¬ 
tity K . The results in the presence of the random strain¬ 
ing are more interesting. We now show that an appropriate 
average of Eq. (I2.76cl) over random b gives a result for short- 
scale energy lossjn accord with Eq. (127741) . Upon recalling 
the definition of b [Eq. (12.621) 1. and noting that the formula 
(I2.76cl) is even in 6, one sees that R(t) = R(0) + 0(t 2 ); 
thus £ > oc R(t) vanishes at t — 0. This is not in con¬ 
flict with formulas like (l2T74ll . however, because those fol¬ 
low from a Markovian closure; one must therefore consider 
times greater than the autocorrelation time r a c of the strain¬ 
ing and evaluate the coarse-grained derivative lim^igi dt £, 
where ‘0’ implies the restriction t r a c- A useful general 
formula for (R(t )) for arbitrary statistics of a (assumed to 
be stationary) seems difficult to obtain; however, one may 
extract the short-time result by expanding 

m = (l+a)- 1 - ^_P(t) + 0(fe 4 ). (2.78) 


One has ( b 2 (t )) = J^dt J^dt' (a(^)a(f , )} « 2(a 2 )r ac t for 
t 3> r ac , which is a standard diffusion law. Thus the coarse¬ 
grained time derivative is 


d(R) 

dt 


i=0 


- 4 ( =| ) (“Vo 


(2.79) 


From £ > = K$ 2 M > and using the fact that J\f > is con¬ 
served, one finds 


d£> _ 2 d{R) 

dt 0 dt 


AT > =- 4 



(a 2 )T ac JV > . 


(2.80) 


To compare this result with Eq. (l2T74ll . we observe that in 
the present model we are assuming that long-wavelength 
straining dominates, so we should assume that 8 k kq is in¬ 
dependent of k. Also, the derivative that is required in 
Eq. (12.74D is explicitly 


d_ 

dk 




(2.81) 


Since the model contains a single K, we take the isotropic 
spectrum A f k = (2n) 2 k~ l 5(k — Kq)M > . One then obtains 
exact agreement between Eqs. (iTSOll and (l2T74ll by replac¬ 
ing r ac by 6 q and choosing 

(a 2 ) = 2 (j^j Af<. (2.82) 


This is nothing but the mean-square strain q 2 U(q); the fac¬ 
tors of q 2 /q 2 correct jV q by removing compressional energy: 
(^Kl 2 ) = 2 q 2 (q 2 /q 2 )£< = 2(q 2 /q 2 ) 2 M<. 


Kraichnan’s model and its generalization assume that 
long-wavelength straining dominates. In general, that is not 
necessarily the case. If short-wavelength decorrelation domi¬ 
nates 9 k q , one must ask whether the factor of R k under 
the k derivative in Eq. (12.741) is all or partly canceled by 

Cf 

the k dependence of 9kkq- For the Galilean-invariant rf k at 
large fc, one can show that in the absence of linear frequen¬ 
cies oc A|(? 7 ^) _1 , or y k oc R k . In the presence of linear 
frequencies, a dependence on R k remains as well, although 
the general case is somewhat complicated. In any case, the 
fact that 8 q kk oc (? 7 fc ) _1 at large k means that the result 
T k oc k k~ l (R k 8 q kk) depends less strongly on R k than 
but is not independent of R k . Clearly the basic conclusion 
that the energy transfer to the large scales is controlled by 
the k dependence of R k 9 q kk still holds. 

2-4-A.3 Relation to zonostrophic instability 

Let us consider the relation between these results and zonos¬ 
trophic instability. For the general anistropic case, if in the 
CE2 zonostrophic instability the A in the denominator of 
Eq. (E3) were replaced by an inverse triad interaction time, 
then Eq. (12.81) in the small q limit agrees with Eq. (12.701) . 
9 does not appear naturally in Eq. (12781) because the CE2 
closure omits eddy damping 77 *.; a more sophisticated closure 
should retain it. A consequence is that the zonostrophic dis¬ 
persion relation derived from CE2 is not correct in all details. 
Nevertheless, we expect that many of its qualitative predic¬ 
tions are robust. The close connection between zonostrophic 
instability and the results derived in this appendix show the 
relevance of the physical mechanism discussed here. In ad¬ 
dition to this physical picture, our discussion has elucidated 
the reason behind the appearance of the factor of R k that 
controls the mathematical behavior of the asymptotic re¬ 
sults. 

Appendix 2.4.B Correlation function 
corresponding to a wave 

We consider in this section the one-time, two-point correla¬ 
tion function corresponding to a wave. First we consider the 
general case of a superposition of waves. Let 

ip'(x, y, t) = 2 ^ c k cos (k x x + k v y - w k f + <p k ). (2.83) 

k 

Then, letting ipi = ip'(xi,yi,t) and V4 = ' l P'( x 2 ,y 2 ,t), we 
have 

=££ 2 c k c k /{ cos [^( k x + k' x )x + ( k x — k' x )x 

k k' 

+ + k y )y + ( k y — k y )y — 2 kk ,] 

+ cos [tj( k x - k' x )x + ( k x + k' x )x 
+ h(ky ~ k y )y + (k v + k y )y — 2 kk ,] }, (2.84) 

where x = x\— X 2 ,x = ^(xi+X 2 ), and 2 kk , = (cn k ±uj k /)f— 
(0 k ± 4>k'). Using a zonal average, the correlation function 
is obtained by integrating over x with x held fixed: 

'&( x ,y\y)=j-[ dx\ x -ip[ Ip 2 , (2.85) 

-for Jo 
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The first cosine vanishes unless k' x = k x , while the second 
cosine vanishes unless k' x = —k x . For simplicity assume all 
the k x ,k x > 0. Then we are left with 

^{x,y | y) = EE 2 c k c k / cos [k x x + ^(k v + k' y )y 

k K 

+ ( k y - k' y )y - (oj k - uw)t + (j > k - 0 k /]. 

( 2 . 86 ) 

If we separate out in the sum the terms for which k' y = k y , 
then we have 


Appendix 2.4.C Dispersion Relation for Arbitrary 
Primary and Arbitrary Secondary Wave 


We show here that for an arbitrary primary wave and arbi¬ 
trary secondary wave, exact agreement is obtained between 
the dispersion relations from CE2 and from the 4-wave mod- 
ulational instability. This generalizes Section [2.4.2.21 which 
shows agreement in the special case where the primary wave 
has p y = 0 and the secondary wave has q x = 0. 

The 4-wave modulational insta bility has a dispersion re¬ 
lation ( Connaughton et al. I teoiolf 12 ! 


y I y) = E 2c k cos {kxX + k y y) 

k 

+ E E 2 c k c k / cos[fc x -a; + \{k y + k' y )y 

k k' y ^k y 

+ (ky - ky)y- (tu k - w k /)t + 0 k - M- 

(2.87) 


id 2 + L d 2 )\-il3q x =Vo|P x q| 2 {P 2 -<1 2 ) 

( 2 2 

E±ZZ _ 

(P+ + L d 2 )(A - iio) - i/3(p x + q x ) 

2 2 \ 

+_ „ p -~ p _) 

( pL + L d “)(A + ioj) + i/3(p x — qx) J 


(2.92) 


It can be verified by substitution that this is a solution to 
the unforced, undamped CE2 equations without zonal flow, 
dtW = 2 Pdydydx’i’ (and using w k = —k x /3/k 2 ). The first 
term of Eq. (1071) . which corresponds to the covariance of 
individual waves, is unchanging in time and homogeneous 
in space. But in the second term, waves with different k y 
give rise to a correlation function that oscillates in time and 
has y dependence. This is a manifestation of the coherent 
beating between waves. 

One can imagine using another averaging procedure in¬ 
stead of the zonal average. With the zonal average, the 
only coherent structures allowed are zonally symmetric. 
One might also want to investigate zonally asymmetric 
structures, which preclude s the use of a zonal average 
( Bakas and Ioannoul . l2013bl f . To study these more general 
coherent structures, the correlation function can be defined 
using a coarse graining in time or space (this approach typ¬ 
ically requires the mean field and fluctuations to obey a 
scale-separation assumption) or an ensemble average. 

To illustrate an alternate derivation for a single wave, let 

V>'(x) = iPo ( e lpx ~ iult + e- ipx+ia,t ) . (2.88) 

Then 

it it 1 2 ( 2ipx —2iu)t , ipx . — ip-x , — 2ipx 2 iujt\ 

W 1 V 2 = Wo ( e e +e^+e+e ^ e I. 

(2.89) 

At this point, a coarse graining in time over an intermediate 
time between cj -1 and the timescale of the coherent struc¬ 
ture eliminates the oscillating terms. Equivalently, one could 
perform a coarse graining in space over an intermediate scale 
between p _1 and the size of the coherent structure. Then, 
one obtains 

T = V>o (e ip x + e" lp x ) . (2.90) 

This T is homogeneous (independent of x) . Its Fourier trans¬ 
form is 

^H{kx,k y ) = (27r) 2 i)>o [<5(k - p) + <5(k + p)]. (2.91) 

The inclusion of the mode at —p as well as the mode at p 
is essential and arises from the reality condition. 


where p± = p ± q and lo = —/3p x /(p 2 + LT 2 ). 

To allow for an arbitrary secondary wave within 
the CE2 formalism , we u se the recent formulation of 
iBakas and Ioannoul 1 2013bl lch. That formulation allows 
for coherent structures of arbitrary spatial dependence 
rather than restricting to zonally symmetric q x = 0 struc¬ 
tures. Their formulation also assumed infinite deformation 
radius, though that could be modified. The dispersion re- 
lation in the small forcing and small dissipation limit is 
( Bakas and Ioannoul . l2013dr^l 


A q 2 - i/3q x = 


where 


f dk x dk y N / 

J ~(2^D V 


1 ~^2 )W H (k x ,k y ), (2.93) 


N — 2 (k x q y — kyq x )^q x q y ( k x + -^-) ~ (ky + ^) 

+ (dy ~ dx) (k x + ^) ( k y + ^-) j, (2.94) 
D = A k 2 k\ - i iq x p [fc 2 + k\ j 

+ 2 if) (k x + - 7 ^) [(foe + -|-) dx + (ky + -|-) dy^ I 


(2.95) 


and k+ = k + q. As in Section 12.4.2.21 the appropriate 
background spectrum to correspond with that of Eq. (127921) 
is Wjj = (27r) 2 i/)gp 4 [<S(k — p) + <5(k + p)]. With sufficient 
algebra, it is possible to show that Eq. (12.931) reduces exactly 
to the L d 2 = 0 limit of Eq. (12.921 . The key is in recognizing 
that 


D = k 2 


N — (k x q y k y q x ) (/i)fj_ k ) , 
^A + - i/3{k x + q x ) 


(2.96) 

(2.97) 


12 This formula correct s a typographical error in Eq. (5.1) of 
I Connaughton et al.l d2010l) . 

13 There is a seeming fac tor of 27 t different from the formula in 
IBakas and Ioannoul (12013d) but that is merely due to the choice 
of Fourier transform convention. 
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2.5 Emergence of non-zonal coherent 
structures — Ioannou & Bakas 
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